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Abstract. We propose a Fast Marching based implementation for com¬ 
puting sub-Riemanninan (SR) geodesics in the roto-translation group 
SE(2), with a metric depending on a cost induced by the image data. 
The key ingredient is a Riemannian approximation of the SR-metric. 
Then, a state of the art Fast Marching solver that is able to deal with 
extreme anisotropies is used to compute a SR-distance map as the so¬ 
lution of a corresponding eikonal equation. Subsequent backtracking on 
the distance map gives the geodesics. To validate the method, we con¬ 
sider the uniform cost case in which exact formulas for SR-geodesics are 
known and we show remarkable accuracy of the numerically computed 
SR-spheres. We also show a dramatic decrease in computational time 
with respect to a previous PDE-based iterative approach. Regarding im¬ 
age analysis applications, we show the potential of considering these data 
adaptive geodesics for a fully automated retinal vessel tree segmentation. 
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1 Introduction 

In this article we study a curve optimization problem in the space of coupled 
positions and orientations M 2 xS' 1 , which we identify with roto-translation group 
SE( 2). We aim to compute the shortest curve y(t) = (x(t),y(t),0(t)) G SE( 2) 
that connects 2 points 7(0) = (xcb2/Cb#o) an d 7(^) = (#i,2/i } 0i) with the re¬ 
striction that the curve is ” lifted” from a planar curve in the sense that the third 
variable 6 is given by Q(t) = arg(x(£) + iy(t )), see Fig. 1. This restriction im¬ 
poses a so-called sub-Riemannian (through the text we denote sub-Riemannian 
as SR) metric that constrains the tangent vectors to the curve to be contained in 
a subspace of the tangent space at every point. This subspace is a plane, which 
differs from point to point, and is the set of all possible tangents to curves in 
SE( 2) that are lifted from smooth planar curves. A SR-metric is a degenerated 
Riemannian metric in which one direction, the one perpendicular to the path 
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in this case, is prohibited (i.e. it has infinite cost). On top of the SR-metric, we 
also consider a smooth external cost, which weights the metric tensor in every 
location and allows for data-adaptivity. 



0 




Fig. 1. Left: the sub-Riemannian problem in SE( 2) can be identified with that of 
a car with two controls (giving gas and steering the wheel). Center: the paths are 
’’lifted” into curves in SE( 2) = R 2 x S' 1 with tangent vectors constrained to the plane 
spanned by the vector fields X\ and X 2 (eq. (pQ)) associated with the controls. Right: 
the SR-spheres (for t — 2, 4 and 6) obtained via the FM-LBR method. 


Essentially, the SR-problem in SE( 2) is that of a car that can go forward, 
backward and rotate (a so-called Reeds-Shepp car) so the possible states of the 
car form a 3D manifold given by the position (x, y) and the orientation 6 of 
the car. Then, admissible trajectories of the car are parametrized by only 2 
control variables associated to the car moving along a straight line (giving gas) 
and to a change of direction (turning the steering wheel). The fact that the car 
cannot step aside infinitesimally imposes the SR-geometry. Finally, the curve 
optimization problem is to find among all possible trajectory between two given 
states, the one with minimal SR-length. 

In image analysis the SR-geodesics in SE( 2) have been proposed in [4] as 
candidates for completion of occluded contours. Here, the geometrical structure 
is used as a model for the functional architecture of the primary visual cortex, 
extending the work in m to the SE( 2) group. This model has proven to be 
valuable in numerous applications (3141516] , and it becomes powerful when com¬ 
bined with the orientation score theory |6] that allows for an invertible stable 
transformation between 2D images and functions on the SE(2) group. The main 
advantage of considering this space of positions and orientations is that inter¬ 
secting curves are automatically disentangled, and therefore the processing in 
the extended domain naturally deals with complex structures such as crossing. 

Sub-Riemannian geodesics in the uniform cost case (the same cost for all 
the SE( 2) elements) were studied by several authors (e.g. (5113] 4. Recently, a 
wavefront based method for computing SR-geodesic that also deals with the 
non-uniform cost case has been proposed in pQ . This method is an extension to 
the SR-case of a classical PDE-approach for computing cost-adaptive geodesics 
used in computer vision, where the metric tensor is induced by the image it- 
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self. The main idea is to consider propagation of equidistant surfaces described 
by level sets of the viscosity solution of an eikonal equation, while subsequent 
backtracking gives the geodesics. In order to solve the eikonal equation, the 
authors rely on a computationally expensive iterative approach based on a left- 
invariant finite-difference discretization of the PDE combined with a suitable 
upwind scheme. 

Here, again in the spirit of computer vision methods, we aim to compute 
the SR-distance map using a Fast Marching method M- This technique, closely 
related to Dijkstra’s method for computing the shortest paths on networks, al¬ 
lows for a significant speed up in the computation of the eikonal equation’s 
viscosity solution. The main difficulty in our case is that classical solvers are 
unable to deal with the extreme (degenerated) anisotropy of the SR-metric. Re¬ 
cently, a modification of the Fast Marching method using lattice basis reduction 
(FM-LBR) that solves this problem was introduced in [9] (code available at 
https://github.com/Mirebeau/ITKFM). Then, the purpose of this paper is to 
show how the SR-curve optimization problem can be numerically solved using 
the FM-LBR method. The key aspects to consider are a Riemannian relaxation of 
the SR-problem and expressing the resulting metric tensor as a matrix-induced 
Riemannian metric in a fixed Cartesian frame. We develop these ideas in the 
Theory section. Then, two experiments are presented. The first one considers 
the uniform cost case (C = 1) and shows that the FM-LBR based method pre¬ 
sented here outperforms the iterative implementation in pQ in terms of CPU time, 
while keeping a similar accuracy. The advantages of considering data-adaptive 
SR-geodesics for extracting blood vessels in retinal images are illustrated in the 
second experiment. 

2 Theory 

Problem formulation. Let g = (x,y,0) be an element of SE( 2 ) = M 2 x5' 1 . 
A natural moving frame of reference in SE( 2 ) is described by the left-invariant 
vector fields {Xi, X 2 , X 3 } spanning the tangent space at each element g: 

X\ = cos Od x + sin Od y , X 2 = do, X% = — sin Od x + cos Od y . (1) 

The tangents j(t) along curves y(t) = (x(t),y(t),6(t)) E SE(2 ) can be written 
as 7 (t) = E?= x i\-y(ty Only the curves with u 3 = 0 are liftings of planar 

curves (see fig. [I]). Then, the tangents to curves that are liftings of planar curves 
are expressed as y(t) = u 1 (t) Xi\^^ +u 2 (t) X 2 \ 1 ^ and they span the so-called 
distribution A = span {Ad, Ad}. Now, the triplet (SE(2), A, Go) defines a sub- 
Riemannian manifold with inner product Go given by: 

Go( 7 , 7 ) = C(j) 2 (/3 2 |x cos sin (9 | 2 + \ 6 \ 2 ) . ( 2 ) 

In view of the car example (see Fig. [Tj) , the parameter /3 > 0 balances between 
the cost of giving gas (the X\ direction) and the turning of the wheel (the X 2 
direction). A smooth function C : SE( 2 ) — [5, 1], 5 > 0 is the given external cost. 
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In order to keep the notation sober during this paper, we do not indicate the 
dependence of Go on the cost C nor that it also depends on the curve 7 . Note that 
the subindex 0 in the metric tensor recalls the SR-structure imposed by allowing 
null displacement in the X 3 direction (i.e. infinite cost for the car to move aside). 
This choice of notation will become clear later the text. Now, optimal paths 7 of 
the car in our extended position orientation space are solutions of the following 
problem: 


W{g) = do{g,e) = min j ^ ||iV)llo dt 

(3) 

where e = (0, 0, 0) is the origin, where || 7 (t)|| 0 = yj Go(7(£), 7(t)) is the SR-norm 
and do is the SR-distance on SE( 2). 

Riemannian approximation. It is possible to obtain a Riemannian ap¬ 
proximation of the SR-problem by expanding the metric tensor in eq. (j2j) to: 

G e ( 7 , 7 ) = £ 0 ( 7 , 7 ) + e _2 C(7) 2 /3 2 |isin6» —ycos6»| 2 , (4) 

where e determines the amount of anisotropy between X% and A. This definition 
bridges the SR-case, obtained at the limit e ^ 0, with the full Riemannian metric 
tensor when e = 1 (isotropic in the spatial directions X\ and X 3 ). Actually, it 
is easy to verify that if C = 1 and (3 = e = 1 , then Gi( 7 , 7 ) = \x \ 2 + \y \ 2 + \0\ 2 . 
Also, by replacing Go with G e in eq. (|3j) we can construct a SE( 2 ) Riemannian 
norm || • || e and a SE( 2) Riemannian distance d e satisfying lim || • || c = || • ||o and 

e^O 

limd e =' do. 

e^O 


7 € A, 7 ( 0 ) = e ,7 (T) = g,T > 0 



f = o.l f -> 0 



Fig. 2. Each ellipsoid represents the Tissot’s indicatrix of the metric G e at different 
elements g G SE(2) (for the case C(g) = 1 and (3=1). The parameter e in eq. Q bridges 
the Riemannian case with the SR-one. When e = 1 each direction has the same cost. 
At the limit e 0, the direction X 3 has infinite cost and the distribution A appears. 


Solution via the eikonal equation. Now we can present the eikonal system 
that solves the problem (j3]) by computing the distance map W(g ) as proved in \ 1 \. 
Following [J], let us introduce some differential operators that will simplify the 
notation of the remaining equations. Let U be a smooth function U : SE( 2 ) M. 
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The Riemannian gradient V e , computed as the inverse of the metric tensor G e 
acting on the derivative, is given by: 

V 6 U := G^dU = C~ 2 r 2 {X 1 U)X 1 + C~ 2 {X 2 U)X 2 + e 2 C“ 2 /3“ 2 (X 3 C)X 3 . (5) 
Then, the norm of the gradient V e is given by: 

||V £ £/|| e = v / C- 2 /3- 2 |X 1 C | 2 + C- 2 |X 2 C | 2 + e 2 C- 2 /3 - 2 |X 3 C| 2 . ( 6 ) 


Thus, the eikonal system that characterizes the propagation of equidistant sur¬ 
faces reads as: 

||V e W(g)|| e = 1, if g^e, 

W(e) = 0. 


( 7 ) 


When e ^ 0 this system becomes the SR-eikonal system in [I] eq.3] where it was 
proved that the unique viscosity solution is indeed the geodesic distance map 
from the origin W(g) = d e (g,e). Then, SR-geodesics are the solutions 7 b(t) of 
the following ODE system for backtracking: 


f i b (t) = t e [o,t] 

\lb(0) = 9- 


(8) 


The metric matrix-representation in the Cartesian frame. A sym¬ 
metric matrix M e representing the anisotropic metric in the frame {ckc, dy , 80} 
can be obtained by a basis transformation from the varying frame {Xi, X 2 , X 3 } 
(see jU Sec. 2 . 6 ]): 


f cos 6 0 — sin 0 N 
M e = | sin 6 0 cos 0 
0 1 0 . 


'C 2 p 2 0 0 

0 C 2 0 
0 0 e~ 2 C 2 /3 2 


f cos 0 0 — sin 6 N 
sin 0 0 cos 0 
0 1 0 


( 9 ) 


Here the diagonal matrix in the middle encodes the anisotropy between the Xi 
directions while the other 2 rotation matrices are the basis transformation. Note 
that the columns are the coordinates of the varying frame in the fixed frame, 
e.g. Xi = cos Odx + sin Ody + 0 80. Then, the metric tensor can be written as 
G e ( 7 , 7 ) = j(t)M e j(t) , with 7 (t) = (x(t),y(t),6(t)). Using this convention, the 
eikonal system © in the fixed frame can be expressed as: 


/ X T W(g)M- 1 XW(g) = 1, if g ± e, 
\W(e) = 0, 


( 10 ) 


where V = (d x , d y ,do) T is the usual Euclidean gradient. The same holds for the 
backtracking equation © which writes: 


r 7 b(t) = -M-'XW^t)), t e [0, T] 
17b(0) = g. 


(ll) 


Note that when approaching the SR-case, the lim M e does not exist but the 

c.R) 

lim is well defined in Eq. (HU). 

e^O 
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Anisotropic Fast Marching. We can now immediately identify Eq. G9 
with 19, eq. 0.1] and then solve the eikonal system via the FM-LBR method. Our 
empirical tests show that e = 0 . 1 , which gives an anisotropy ratio k = 0.01 (see 
[9j eq. 0.5]), is already a good approximation of the SR-case and is the value 
used in the following experiments. 

3 Experiments and Applications 

Validation via comparison in uniform cost case. The exact solutions of the 
SR-geodesic problem for the case C = 1 are known (see m for optimal synthesis 
of the problem). Therefore, and similar to what is done in [I], we consider this 
case as our golden standard. Here, we want to compare both the computational 
time and the accuracy achieved in the calculation of the discrete SR-distance 
map W(g ), which is the solution of the eikonal system © when e ^ 0 . 

Let us set /3 = 1 and consider a grid Q s = {(xi, yj, 0 k) G SE(2)} with uniform 
step size s, where xi = is , yj = js, 0k = fcs, the indices i,j,k G N such that 
\xi\ < 27r, \yj | < 27r and — tt+s < 0k < 7r. Then we compute the discrete geodesic 
distance map W(g) on Q s using the iterative method in [T] and the FM-LBR. In 
order to measure the accuracy of the achieved solutions we follow the method 
explained in detail in [T. There, by solving the initial value problem from the 
origin e, a set of arc length parametrized SR-geodesics is computed such that 
SR-spheres of radius t are densely sampled. Then, the endpoints g = (x, y , 0) of 
each geodesic is stored in a list together with its length t. Finally, we define the 
max relative error as E^^t) = max(| W(g) —t\/t) where the max is taken over all 
the endpoints g and where the value of W(g) is obtained by bi-linearly interpo¬ 
lating the numerical solutions of eq. m computed on the grid Q s . The results and 
comparisons are presented in Fig. [3j Here we solved the eikonal equation in in¬ 
creasingly finer grids Q s obtained by setting the step size s = 7r/n,n G N + . Note 
that the size of Q s is then (2n + l)(2n + l)(n —1). The graph in Fig. [3^ left) shows 
the comparison of the accuracy achieved in the computation of the SR-sphere of 
radius t = 4 when n increases. The behaviour for SR-spheres of different radius 
is similar. The CPU time is compared in Fig. [^center). The 3rd plot illustrates 
the method for computing the accuracy. The orange surface is the SR-sphere of 
radius t = 4 computed with the FM-LBR method on a grid corresponding to 
n = 101. The points are the geodesic endpoints, their color is proportional to the 
error of the FM-LBR (blue-low, green-medium, red-high error). The first obser¬ 
vation is that even though the iterative method is more accurate, both methods 
seem to have the same order of convergence (the slope in the log-log graphs) 
when the grid resolution increases. This seems reasonable as both methods use 
first order approximations of the derivatives. Also, we hypothesise that the off¬ 
set in favour of the iterative method is due to the Riemannian approximation 
of the SR-metric (i.e. selecting e = 0.1), but this needs further investigation. 
The second key observation is that the CPU time increases dramatically with 
n for the iterative method. Therefore, it is clear that we can achieve the same 
accuracy using the FM-LBR but with much less computational effort, which is 
of vital importance in the subsequent application. 
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Fig. 3. Validation via comparison in the uniform cost case. The experiment (illustrated 
on the left, see the text) shows that even though the iterative method in [l] is more 
accurate we can still achieve with the FM-LBR method better results and with less 
CPU effort, just by increasing the grid sampling. 


Application to retinal vessel tree extraction. The analysis of the blood 
vessels in images of the retina provides with early biomarkers of diseases such 
as diabetes, glaucoma or hypertensive retinopathy [7]. For these studies, it is 
important to track the structural vessel tree, a difficult task especially because 
of the crossovers and bifurcations of the vessels. Some existing techniques man 
rely in considering an extended (orientation and/or scale) domain to deal with 
this issue. Moreover, in pQ promising results were obtained by formulating the 
vessel extraction as a SE(2 ) SR-curve optimization problem with external cost 
obtained through some wavelet-like transformation of the 2D images. In the 
previous experiment, we have shown that our proposed FM-LBR based imple¬ 
mentation computes in practice the same geodesics as the iterative method in \T\. 
Therefore, by simply replacing the eikonal equation solver in [L we can obtain 
the same results but with a dramatic decrease of the computational demands 
(both CPU time and memory). The example in Fig. 0] shows the advantages of 
considering the SR-metric in performing the vessel tree extraction. In this case 
(a patch of size 200x200, with 64 orientations considered) the iterative method 
computed the distance map in approximately 1 hour while the FM-LBR did the 
same in 20 seconds. For the experiments details we refer to pQ, for more examples 
see www.bmia.bmt.tue.nl/people/RDuits/Bekkersexp.zip. 


4 Conclusions 

Over the last decade, some authors mm have shown the advantages of consider¬ 
ing the roto-translation group embedded with a SR-geometry as a powerful, rich 
tool in some image analysis related problems or for the geometrical modelling 
of the visual perception. In our opinion, 2 obstacles have prevented this frame¬ 
work to become more popular amongst engineers: the expensive computational 
demands involved (resulting of considering the extended orientation space) and 
the lack of efficient numerical methods able to deal with the extreme (degener¬ 
ated) anisotropy imposed by the SR-metric. These obstacles are addressed by 
the main contribution of this work, which is solving (up to our knowledge for the 
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Fig. 4. Tracking of blood vessels in retinal images via cost adaptive SR-geodesics (see 
experiment details in p]). Left: the cost obtained from the image (orange indicates 
locations with low cost). Center: Tracking in the full (e = 1) Riemannian case. Right: 
Tracking with the approximated (e = 0.1) SR-geodesics. 


first time) the SR-geodesic problem using a Fast Marching based implementa¬ 
tion. To be able to achieve this, we rely on the FM-LBR solver recently introduce 
in [9] and show that even when relaxing the SR-restriction by a Riemannian ap¬ 
proximation of the metric we achieve excellent numerical convergence, but much 
faster than with the approach in pQ. Regarding the retinal imaging application 
our promising preliminary studies suggest that it is at least feasible to aim for 
a full vessel tree segmentation as the solution of a single optimization prob¬ 
lem, but this requires further investigation. Future work will pursue extension 
of this method to the 3D-rototranslation group SE( 3) and the applications in 
neuroimaging and 3D-vessel segmentation. 
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